#### Calibration error estimatation

## Setup environment

setwd("~/research/coffee-dataverse/src")

do.origspec <- T
do.configonly <- T
source("load.R", encoding="iso-8859-1")

library(rstan)
library(dplyr)

## Collect information form two calibrations

metas <- data.frame()
munis <- data.frame()
for (output.suffix in c("-eq24-orig", "-eq24-orig2")) {
    ## Load data
    param <- get.fit.params(output.suffix)
    load(paste0("../data/combined-betagamma", output.suffix, ".RData"))

    sumstats <- read.csv(paste0("../data/sumstats", output.suffix, ".csv"))
    names(sumstats)[1] <- "fullparam"
    sumstats$param <- gsub("\\[\\d+\\]", "", gsub("_coeff\\[(\\d+)\\]", "_coeff\\1", sumstats$fullparam))
    names(sumstats)[names(sumstats) == 'n_eff'] <- 'neff'

    muniinfo <- read.csv(paste0("../data/combined-muniinfo", output.suffix, ".csv"))

    for (pp in 1:length(param)) {
        print(param[pp])
        if (param[pp] == 'pricememory') {
            name <- 'autoreg_price'
            values <- 1 - la$gamma[, 1, pp]
        } else if (param[pp] == 'autoreg') {
            name <- 'autoreg_planting'
            values <- la$gamma[, 1, pp]
        } else {
            name <- param[pp]
            values <- la$gamma[, 1, pp]
        }
        betas <- la$beta[,, pp]

        mu <- mean(values)
        se <- sd(values)
        lo95 <- quantile(values, .025)
        hi95 <- quantile(values, .975)
        metas <- rbind(metas, data.frame(model=output.suffix, param=name, mu, se, lo95, hi95))

        myparam <- param[pp]
        muniinfo2 <- data.frame(region=muniinfo[, 2]) %>% left_join(subset(sumstats, param == myparam) %>% group_by(region) %>% summarize(neff=mean(neff, na.rm=T), Rhat=mean(Rhat, na.rm=T)))
        muniinfo2$mu <- colMeans(betas)
        muniinfo2$se <- apply(betas, 2, sd)

        munis <- rbind(munis, cbind(model=output.suffix, param=name, muniinfo2))
    }
}

## Ensure that parameters match up
all(metas$param[metas$model == '-eq24-orig'] == metas$param[metas$model == '-eq24-orig2'])

## Prepare labels

labels <- list('death_coeff1'="Die-off EDD effect", 'autoreg_planting'="Planting autoreg.", 'autoreg_price'="Nerlove autoreg.", 'cost_harvest'="Cost of harvest", 'cost_to_planting'="Planting price resp.", 'scale_trend'="Exogenous trend", 'yield_sigma'="Yield std. dev.", 'yield_uncertainty'="Range of yields", 'bearing_sigma'="Bearing std. dev.", 'harvest_sigma'="Harvest std. dev.", 'production_sigma'="Production std. dev.")
symbols <- list('death_coeff1'='\\delta', 'autoreg_planting'='\\rho_a', 'autoreg_price'='\\rho_p', 'cost_harvest'='c', 'cost_to_planting'='\\phi', 'scale_trend'='\\theta', 'yield_sigma'='\\sigma_y', 'yield_uncertainty'='\\Delta', 'bearing_sigma'='\\sigma_b', 'harvest_sigma'='\\sigma_h', 'production_sigma'='\\sigma_q')
for (kk in 1:length(weatherpreds)) {
    labels[[paste0("yield_coeff", kk)]] <- weatherlabels[kk]
    symbols[[paste0("yield_coeff", kk)]] <- weathersymbols[kk]
}

pmetas <- data.frame(Parameter=sapply(metas$param[metas$model == '-eq24-orig'], function(p) labels[[p]]),
                     Symbol=paste0("$", sapply(metas$param[metas$model == '-eq24-orig'], function(p) symbols[[p]]), "$"),
                     Mean1=metas$mu[metas$model == '-eq24-orig'], SE1=metas$se[metas$model == '-eq24-orig'],
                     Mean2=metas$mu[metas$model == '-eq24-orig2'], SE2=metas$se[metas$model == '-eq24-orig2'])

## Print table

library(xtable)

print(xtable(pmetas), include.rownames=F, sanitize.text.function=function(x) x)

## Determine relationship between error and N_eff and Rhat

library(nnls)

pmunis <- data.frame()
for (pp in unique(munis$param)) {
    orig <- subset(munis, param == pp & model == '-eq24-orig')
    addl <- subset(munis, param == pp & model == '-eq24-orig2')
    combo <- orig %>% left_join(addl, by='region', suffix=c('', '2'))

    combo$logsqrdiff <- log((combo$mu - combo$mu2)^2)

    mod <- lm(logsqrdiff ~ log(neff) + log(Rhat), data=combo)

    pmunis <- rbind(pmunis, data.frame(param=pp, const=mod$coeff[1], logneff=mod$coeff[2], logneff.se=sqrt(vcov(mod)[2, 2]),
                                       logRhat=mod$coeff[3], logRhat.se=sqrt(vcov(mod)[3, 3])))
}

get.stars <- function(mu, se) {
    ifelse(mu/se > 3.29, "***", ifelse(mu/se > 2.58, "**", ifelse(mu/se > 1.96, "*", ifelse(mu/se > 1.65, '.', ''))))
}

## Print table

pmunis2 <- data.frame(Parameter=sapply(metas$param[metas$model == '-eq24-orig'], function(p) labels[[p]]),
                      Symbol=paste0("$", sapply(metas$param[metas$model == '-eq24-orig'], function(p) symbols[[p]]), "$"),
                      Constant=format(pmunis$const, digits=3),
                      `log(Neff)`=paste0(format(pmunis$logneff, digits=1, scientific=F), "$^{", get.stars(pmunis$logneff, pmunis$logneff.se), "}$"),
                      `(SE)`=paste0("(", format(pmunis$logneff.se, digits=1, scientific=F), ")"),
                      `log(Rhat)`=paste0(format(pmunis$logRhat, digits=1, scientific=F), "$^{", get.stars(pmunis$logRhat, pmunis$logRhat.se), "}$"),
                      `(SE)`=paste0("(", format(pmunis$logRhat.se, digits=1, scientific=F), ")"))

print(xtable(pmunis2), include.rownames=F, sanitize.text.function=function(x) x)
